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When liquid-crystal elastomers are prepared without any alignment, disordered polydomain struc¬ 
tures emerge as the materials are cooled into the nematic phase. These polydomain structures have 
been attributed to quenched disorder in the cross-linked polymer network. As an alternative expla¬ 
nation, we develop a theory for the dynamics of the isotropic-nematic transition in liquid-crystal 
elastomers, and show that the dynamics can induce a polydomain structure with a characteristic 
length scale, through a mechanism analogous to the Cahn-Hilliard equation for phase separation. 


Liquid-crystal elastomers are remarkable materials 
that combine the elastic properties of cross-linked poly¬ 
mer networks with the anisotropy of liquid crystals !fj. 
Any distortion of the polymer network affects the orien¬ 
tational order of the liquid crystal, and any change in the 
magnitude or direction of liquid-crystal order influences 
the shape of the polymer network. Hence, these elas¬ 
tomers are useful for applications as actuators or shape¬ 
changing materials. 

For many applications, it is necessary to prepare mon¬ 
odomain liquid-crystal elastomers. In practice, this can 
be done by applying a mechanical load or other align¬ 
ing field while crosslinking [2]. Surprisingly, elastomers 
prepared without an aligning field do not form mon¬ 
odomains. Rather, they form polydomain structures 
with nematic order in local regions, which are rnacro- 
scopically disordered. These polydomain structures have 
been seen in many experiments, using a wide range of 
techniques EH3- Indeed, a recent polarized light scat¬ 
tering study shows that liquid-crystal elastomers evolve 
toward a state of increasing disorder as the isotropic- 
nematic transition proceeds, unless the disorder is sup¬ 
pressed by a gradually increasing load 0. 

One important issue in the theory of liquid-crystal 
elastomers is how to understand the polydomain state. 
Several theoretical studies have attributed this state to 
quenched disorder in the polymer network, which can be 
understood by analogy with spin glass theory GUI!- Ef¬ 
fects of quenched disorder have further been modeled and 
visualized through numerical simulations HHM]. More 
macroscopic theories have shown that the resulting poly¬ 
domain structure has profound consequences for the ma¬ 
terial’s effective elasticity [221123] ■ 

The purpose of this paper is to suggest a different 
mechanism for the origin of the polydomain state, not re¬ 
lated to quenched disorder. We develop a theory for the 
dynamics of the isotropic-nematic transition in liquid- 
crystal elastomers, in which growing nematic order is 
coupled to elastic strain. We explore this theory in two 
dimensions (2D), using two models for dynamic evolu¬ 
tion of nematic order and strain. The theory shows that 
dynamics can itself select a characteristic length scale 
for a disordered polydomain structure, through a mech¬ 
anism similar to the Cahn-Hilliard equation for phase 


separation. In particular, the theory predicts formation 
of structures with the form shown in Fig. |T| We sug¬ 
gest that this mechanism may play a role in formation 
of polydomain liquid-crystal elastomers, in addition to 
quenched disorder. 

In the theory of 2D liquid-crystal elastomers, nematic 
order is described by the tensor order parameter Q a p{r), 
and elastic distortion of the material by the displacement 
vector u(r). In terms of displacement, the strain tensor 
is defined as e a p = \[d a up + dpu a + (d a u 1 )(d/3U 1 )]. The 
free energy can be expressed in terms of Q a p and e a p as 


F= d 2 r 
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Here, the first two terms are the Landau-de Gennes ex¬ 
pansion for the free energy in powers of the order tensor. 
The coefficient a = a'(T — T 0 ) is assumed to vary linearly 
with temperature, while & is a positive constant. The 
third term is the Frank free energy for spatial variations 
in the order tensor, assuming a single Frank coefficient L. 
The fourth and fifth terms are the elastic free energy in 
terms of the strain tensor, with Lame coefficients A and 
/r. The final term is the coupling between nematic order 
and strain, with coefficient F. 

If there were no coupling between nematic order and 
strain, F = 0, the system would have an isotropic- 
nematic transition at a = 0, corresponding to temper¬ 
ature T 0 . With coupling F/0, the transition is shifted 
upward to a = V 2 / (2 fi), corresponding to the higher tem¬ 
perature Tin = T 0 +V 2 /(2fj,a'). Above that temperature, 
the state of minimum free energy is uniformly isotropic, 
with Q a p = 0 and e a p = 0. Below that temperature, 
at a = V 2 /{2fi) — 5a, the state of minimum free energy 
becomes uniformly nematic, with alignment along a ran¬ 
domly selected director h. In this state, the order ten¬ 
sor is Q a p = S(2n a np — 5 a p), where the magnitude of 
nematic order is S = \J5a/(2b). This state extends uni¬ 
formly along the director, with strain e a p = [F/(2 ft)]Q a p- 

Now suppose we begin in the isotropic phase, and 
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FIG. 1. Visualization of liquid-crystal elastomer structures 
calculated here. The orientation and eccentricity of ellipses 
represents local nematic order, (a) Single wave in nematic 
order and displacement, (b) Superposition of two perpen¬ 
dicular waves, forming a square lattice, (c) Superposition of 
three waves with random directions, amplitudes, and phases. 


rapidly cool to a temperature slightly below Tin- At this 
low temperature, nematic order and strain both begin to 
grow dynamically. We ask: Does the dynamic process 
lead to the state of minimum free energy, with uniform 
Q a p and e a ^? Alternatively, does it lead to a different, 
nonuniform state? 

To answer this question, we develop a model for the 
dynamics of the phase transition. We actually consider 
two models, first simple linear drag and then more real¬ 
istic viscous flow. In both models, we describe four cou¬ 
pled degrees of freedom: Q xx (r,t), Q xy (r,t), u x (r,t), 
and u y (r,t). The remaining components of Q a p( r ,t) 
are fixed because it is a symmetric, traceless tensor, and 
e a p(r,1) can be derived from u(r, t). We cannot take the 


strain tensor components as our fundamental degrees of 
freedom because they are constrained by elastic compat¬ 
ibility; they must all be derivable from the same u(r,t). 

In the simplest model of overdamped dynamics with 
linear drag, the rate of change for each degree of freedom 
is linearly proportional to the force acting on it. Hence, 
the equations of motion are 

dQ xx _ SF dQ xy _ SF 
dt Q SQ XX ’ dt Q SQ xy 

du x SF du v ^ SF 

~dt = “W "i dt = ~ Tu dUy ’ ^ 


where Tq and r u are mobility coefficients. To calculate 
the forces on the right side of those equations, we sub¬ 
stitute the definition of the strain tensor into the free 
energy and take functional derivatives with respect 
to Q a p and u a . We then linearize the equations, assum¬ 
ing that Q a p and u a are both small in early stages of 
nematic ordering. The equations then become 

= r Q [-2 aQ xx + 2LS7 2 Q XX + V(d x u x - d y u y )] , 
dO 

— = r<g [— 2aQ xy + 2U\7~Q xy + V(d x u y + d y u x )^ , 
= r u [(A + n)d x V ■ u + [iV 2 u x - V( V • Q) x \ , 

= T u [(A + fi)d y V ■ u + ti\/ 2 u y - V(V • Q)y] . 

(3) 


To simplify this system of equations, we Fourier trans¬ 
form from position r to wavevector k , then write the 
equations in the matrix form 



^Qxx{k,tF 


(Qxx{k,t)\ 

d 

Ft 

Qxy(h, t) 

u x (k,t) 

= -M(fe) 

Qxy(k, f) 

u x (k, t) 


\ Uy{k,t) y 


\ u y{k-i t) ) 


where M(k) is a 4 x 4 matrix. This matrix equation re¬ 
sembles the Cahn-Hilliard equation for phase separation 
of a binary fluid. At each fc, the matrix M(k) has four 
eigenmodes i, which either grow or decay exponentially 
as e _A d fc ) t , where Aj(fe) is the corresponding eigenvalue 
oiM(k). Note that Ai(fc) < 0 corresponds to exponential 
growth, while Aj(fe) > 0 corresponds to exponential de¬ 
cay. We must determine what grows most rapidly: which 
eigenmode at which wavevector? 

To identify the fastest-growing mode, we choose coor¬ 
dinates such that k is along the x-axis. The matrix then 
simplifies to 


Af(fc) = (5) 

/2Tq(ci + Lk 2 ) 0 -iT Q Vk 0 \ 

0 2T Q {a + Lk 2 ) 0 -iV Q Vk 
iT u Vk 0 T u (A + 2 n)k 2 0 

\ 0 iT u Vk 0 T u nk 2 / 
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FIG. 2. Sample plot of the eigenvalues A±(fc) in the linear 
drag model of dynamics, with parameters a = 0.1 and L = 
/i = V = r q = F u = 1. The largest negative eigenvalue 
corresponds to the fastest-growing mode, which occurs at a 
dynamically selected wavevector. 


We now take the limit of an incompressible material, with 
A —¥ oo. In this limit, u x has a high energy cost, so that 
it decays rapidly, and hence we eliminate it from consid¬ 
eration. In that case, Q xx is not coupled to any other 
degrees of freedom, so it is an eigenmode by itself, with 
eigenvalue 2Fg(a + Lk 2 ). If the system is at a temper¬ 
ature slightly below the isotropic-nematic transition, we 
must have 0 < a < V 2 /{2n). In that temperature range, 
this eigenvalue is positive, so that Q xx decays exponen¬ 
tially. Hence, we also eliminate it from consideration in 
the search for the fastest-growing mode. 

The remaining two modes are linear combinations of 
Q xy and u y , with eigenvalues 

A±(k) =V Q (a +Lk 2 ) + ^T u fik 2 (6) 

± ^[T Q (a + Lk 2 )~- \T u nk 2 } 2 + TqT u VW 

Figure [2] shows a sample plot of these two eigenval¬ 
ues as functions of k. The eigenvalue A + (fc) begins 
at 2Tga when k = 0, then increases with increasing 
k. For temperatures just below the isotropic-nematic 
transition, with 0 < a < V 2 /(2n), it is always posi¬ 
tive and hence represents a decaying mode. By con¬ 
trast, A -ik) begins at 0 when k = 0, then decreases 
into negative values over the range 0 < k < yJ~6aj~L, 
where 6a = V 2 /{2n) — a , and eventually returns to 
positive values for larger k. Over the range in which 
it is negative, it represents an exponentially growing 
mode. To find the fastest-growing wavevector, we min¬ 
imize A_(fc) over k. For temperatures just below the 
isotropic-nematic transition, for small 6a, this wavevec¬ 
tor is fcf as test ~ \j6a/ (2 L), and the corresponding growth 
rate is |A_(/cf as test)| ~ V u \i 2 5a 2 /{2LV 2 ). 

We emphasize that this wavevector is selected through 
a dynamic mechanism. It is not the minimum of the free 
energy (which is a state of uniform nematic order and 


strain). Moreover, it only occurs because of the coupling 
V between nematic order and strain in a liquid-crystal 
elastomer. If these variables were uncoupled (V = 0), 
the matrix M would be diagonal, the isotropic-nematic 
transition would occur at a = 0, and the fastest-growing 
mode below that transition would be k = 0. 

To characterize the fastest-growing mode, we calculate 
the eigenvector of M corresponding to eigenvalue A_ at 
wavevector fcf astes t- This eigenvector represents waves in 
both Q xy and u y (with our assumption that the wavevec¬ 
tor is in the x-direction), and these waves are 90° out of 
phase. Figure [l|a) shows a sample visualization of the 
structure with a single Fourier mode. It has alternat¬ 
ing stripes with the director oriented at ±45° from the 
wavevector, accompanied by displacement perpendicular 
to the wavevector. 

In general, a liquid-crystal elastomer will not have only 
one Fourier mode. Rather, it can include modes with 
wavevectors of magnitude fcf as test hr multiple directions. 
To find a mode in an arbitrary direction, we rotate the 
wavevector, and make a corresponding rotation of Q a p 
and u. We then add up the Fourier modes to find the 
structure. Figure [Tj(b) shows an example with two per¬ 
pendicular waves of equal amplitude, leading to a square 
lattice in the nematic order and the displacement. Fig¬ 
ure [ljc) shows a more realistic example with a superposi¬ 
tion of three waves with random directions, amplitudes, 
and phases. 

The structures in Fig. [l] are similar to structures com¬ 
monly observed in experiments and simulations on active 
nematic liquid crystals [24] . This similarity is reason¬ 
able, because both systems are controlled by couplings 
between orientational order and extension of the mate¬ 
rial. 

The growth of nematic order in a liquid-crystal elas¬ 
tomer can be described by the dynamic correlation func¬ 
tion 


C(\ r — r'\ : t) = (cos 2 [0(r)-0(r')]) t 
— (Qxxij* i tyQxxiy 1QxyiX i ^)Qxy{j* ? ^)) 


= £■ 


} ik-(r—r') 


{\Qxx{k,t)\ 2 + \Q X y(k,t)\ 2 ) . (7) 


This sum is dominated by the fastest-growing mode at 
wavevectors with magnitude fcf as t es tj and hence 

C'(\r — r'\,t) (X e ifc '( T ’- T ’') e 2 l A -( fe fastest)|t 

|^| = ^fastest 

OC Jo(fcfastest> - r'|)e 2|A - (fcf “ test)|t . (8) 


Thus, in the early stages of growth, the correlation func¬ 
tion has the form of Bessel function Jo (^fastest k — ^l), 
with an exponentially increasing magnitude. In later 
stages of growth, the approximation of small nematic 
order ceases to apply, and other types of modeling are 
needed. Even so, the length scale of l/fcf as t es t is estab¬ 
lished from the early stages. 
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The dynamic model presented above has a limitation: 
It assumes that both Q a p(r,t) and u(r,t) have over¬ 
damped dynamics, with drag forces linearly proportional 
to the rate of change of these quantities. This assump¬ 
tion is appropriate for dynamics on a substrate, where 
the dissipation is caused by drag against the substrate. 
However, if there is no substrate, it is reasonable to gen¬ 
eralize the dynamics in two ways: by considering inertia 
for the displacement and by considering viscous dissipa¬ 
tion rather than drag against a substrate. 

For this generalization, we use the equations of motion 


d 2 u a SD SF 

^ dt 2 5u a 8u a ’ 


0 = — 


SD 

SQa0 


SF 

SQafl 


(9) 


Here p is the mass density, which gives inertia for u ; there 
is no inertia for Q a p. Also, D is the Rayleigh dissipation 
function, which can be written as 


D= d 2 


pA a pA a p T ^ 71 D a pB a p -f- r )2A a pB OL p 


( 10 ) 


in terms of the two modes that dissipate energy: A a p = 
\{d a up + dpu a ) is the rate of shear flow, and B a p = 
Qa@ — UziesaQsp + f-spQsa) is the rotation rate of ne¬ 
matic order relative to rotational flow of the material, 
given by lo z = ^e^dyii^. In these expressions, p is the 
viscosity, 71 is the rotational viscosity, and 72 is a dissi¬ 
pative coupling coefficient. 

We combine these expressions to derive the coupled 
equations of motion for Q xx , Q xy , u x , and u y , and lin¬ 
earize the equations assuming these variables are small in 
the early stages of nematic ordering. We then follow the 
same steps as in the previous calculation: Fourier trans¬ 
form from r to k, choose coordinates such that k is along 
the x-axis, eliminate u x by the constraint of incompress¬ 
ibility, and eliminate Q xx because it is an independent, 
exponentially decaying mode. We are left with a matrix 
equation for Q xy (k,t ) and u y (k,t), 




( 4 7 i i^ 2 k\ (Q xy \ , x 

\-i72fc 77 k 2 ) \ u y ) 

(2(a + Lk 2 ) -iVk\fQ xy \ 
iVk pk 2 ) \ u y J ' 


Next we Fourier transform from time t to frequency w, 
and obtain 


/2(a + Lk 2 ) — 4foia; —iVk + 7 2 wfc \ fQ xy \ 
iVk — ^uik pk 2 — irjuik 2 — pui 2 J \u y J 

(! 2 ) 

In this matrix equation, there are two couplings be¬ 
tween Q xy and u y : the elastic coupling V and the dis¬ 
sipative coupling 72 . For simplicity, we set 72 = 0 and 
consider only the elastic coupling. 


-Im[oj(k)l 



FIG. 3. Sample plot of the mode structure in the generalized 
model of dynamics, with inertia and viscosity. Parameters 
are a = 0.1, L = p = V = r/ = 1, 71 = 0.25, and p = 20. The 
quantity — Im[u;(fcy| is the exponential decay rate, equivalent 
to A±(fc) in Fig. [ 2 ] The largest positive value of Im[aj(fc)] 
corresponds to the fastest-growing mode. 


The matrix equation only allows nontrivial Q xy and u y 
if the determinant of the matrix is zero. Hence, we set the 
determinant to zero and solve for the allowed frequencies 
w. Because the determinant is a cubic function of w, there 
are three solutions. Expanding to first order in 1/p, the 
solutions are 


i(a + Lk 2 ) 

+ 


Z7ifo 2 fc 2 
p(a + Lk 2 ) 2 ’ 

n\ ^ ik2 

“W = ± ‘v - MPPm - 7? 


(13) 


V + 


'yiV 2 


(■a + Lk 2 ) 2 


Here, a real part of oj represents oscillation, a negative 
imaginary part represents exponential decay, and a pos¬ 
itive imaginary part represents exponential growth. 

The solution wo(fc) is a purely damped mode. Whether 
the system is in the isotropic phase, a > V 2 /(2p), or 
slightly in the nematic phase, 0 < a < V 2 /{2p), this 
mode decays exponentially. 

The modes uj± ( k ) depend on whether system is in 
the isotropic or nematic phase. In the isotropic phase, 
a > V 2 /(2p), these modes are damped sound waves, with 
both oscillation and exponential decay. By comparison, 
when the system is cooled slightly into the nematic phase, 
0 < a < V 2 /(2p), these modes change into pure expo¬ 
nential growth or decay. One of the modes has a negative 
imaginary part for all k , corresponding to decay, but the 
other mode has a positive imaginary part for a range of 
k , corresponding to growth. 

Figure [3] shows a sample plot of the mode structure 
in the nematic phase. We can see that the ujo(k) and 
u)-(k ) modes are decaying for all fc, but the w+(fc) mode 
is growing for a range of k. In this respect, it resem¬ 
bles the growing mode in the linear drag model of dy¬ 
namics, shown in Fig. [2] In the limit of high p, the 
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range of exponential growth is 0 < k < \J8a/L, and 
the fastest-growing wavevector is /cf aste st ~ \JSa/ (2L), 
where Sa = V 2 / (2 fi) — a. These results are equivalent to 
corresponding results for the linear drag model. 

Hence, the generalized model of dynamics (with inertia 
and viscosity) leads to the same conclusion as the linear 
drag model: The dynamic mechanism of the isotropic- 
nematic transition selects a fastest-growing wavevector. 
This fastest-growing wavevector is not the minimum of 
the free energy, and it only occurs because of the cou¬ 
pling between nematic order and strain. We expect mod¬ 
ulations with this wavevector to grow in liquid-crystal 
elastomers cooled below the isotropic-nematic transition, 
leading to structures with the form shown in Fig. [l] 

To be sure, both models of dynamics presented here 
apply only to early stages of growth of nematic order. 
In later stages, as nematic order becomes more estab¬ 
lished, we cannot assume that Q a p{ r ^) and u(r,t) are 
small. In that case, our linearization of the equations of 
motion breaks down, and the dynamics must be stud¬ 
ied through other techniques, such as numerical simula¬ 
tion. Hence, we cannot be sure whether the polydomain 
structure will persist into longer time, or will eventually 
coarsen into a uniform structure. In those stages of dy¬ 
namics, quenched disorder may lock in the polydomain 
structure at the length scale given by dynamics, and pre¬ 
vent it from coarsening away. 

In conclusion, we have shown that dynamic evolution 
of nematic order can induce a polydomain state with 
a characteristic length scale, in the early stages of the 
isotropic-nematic transition. This mechanism should be 
considered, along with quenched disorder, in studies of 
polydomain liquid-crystal elastomers. 

We thank D. R. Nelson and M. Y. Pevnyi for helpful 
discussions. This work was supported by NSF Grant 
DMR-1409658. 
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